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ABSTRACT 

In a previous study, we formulated a framework of the entropy -based equilibrium statistical 
mechanics for self-gravitating systems. This theory is based on the Boltzmann-Gibbs entropy 
and includes the generalized virial equations as additional constraints. With the truncated 
distribution function to the lowest order, we derived a set of second-order equations for the 
equilibrium states of the system. In this work, the numerical solutions of these equations are 
investigated. It is found that there are three types of solutions for these equations. Both the 
isothermal and divergent solutions are thermally unstable and have unconfined density pro- 
files with infinite mass, energy and spatial extent. The convergent solutions, however, seem 
to be reasonable. Although the results cannot match the simulation data well, because of the 
truncations of the distribution function and its moment equations, these lowest-order conver- 
gent solutions show that the density profiles of the system are confined, the velocity disper- 
sions are variable functions of the radius, and the velocity distributions are also anisotropic in 
different directions. The convergent solutions also indicate that the statistical equilibrium of 
self-gravitating systems is by no means the thermodynamic equilibrium. These solutions are 
just the lowest-order approximation, but they have already manifested the qualitative success 
of our theory. We expect that higher-order solutions of our statistical-mechanical theory will 
give much better agreement with the simulation results concerning dark matter haloes. 

Key words: methods: analytical - galaxies: kinematics and dynamics - cosmology: theory - 
dark matter - large-scale structure of Universe. 



1 INTRODUCTION 

This is a follow-up work to our series of investigations of the 
statistical mechanics of self-gravitating systems. It has long been 
realized that the conventional methods of statistical mechanics 
of short-range interaction systems cannot be directly applied to 
the study long-range self-gravitating systems. Hence, it is neces- 
sary to return to the fundamentals of statistical mechanics and de- 
vel op special technique s to handle the l ong-range nature of grav- 
itv jPadmanabhanll 19901: IChavanisll2q06l) . In our earlier two works 
( iHe & Kandl20ia : lKang & Hdl201ll ). we performed preliminary 
investigations into the statistical mechanics of coUisionless self- 
gravitating systems, and found that both the concept of entropy 
and the entropy principle are still valid for self-gravitating systems, 
which implies that an entropy-based statistical mechanics of self- 
gravitating systems is feasible. 

In|He & Kane ( 201 1), we also investigated the stability of the 
statistical equilibrium of spherically symmetric coUisionless self- 
gravitating systems. By calculating the second variation of the en- 
tropy, we found that perturbations of the relevant physical quan- 
tities should be classified as long- and short-range perturbations. 
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which correspond to the long- and short-range relaxation mecha- 
nisms, respectively. We showed that the statistical equilibrium state 
with the most probable distribution of self-gravitating systems is 
neither maximum nor minimum, but is complex saddle-point en- 
tropy state, and hence differs greatly from the case of an ide al gas. 
Violent relaxation should be divided into two phases (see fSokeJ 
1 1996) . The first phase is the entropy-production phase, while the 
second phase is the entropy-decre asing phase. The saddle-point en- 
tropy state is consistent with both Antonov ( 1962)'s proof and Bin- 
ney's argument (see Binney & Tremaine 2008). 

Based on these findings, iHel l l2012h formulated a system- 
atic theoretical framework of the statistical mechanics for spher- 
ically symmetric coUisionless self-gravitating systems. It was 
demonstrated that the ergodicity is invalid globally for the whole 
self-gravitating system but can be re-established locally, and the 
Boltzmann-Gibbs entropy was re-derived, which should be the 
correct entropy form for coUisionless gravitating systems. This 
framework i s significantly diffe rent from previo u s works (e.g . 
Lvnde n-Bell 1967; Shu 1978: Kull et alj 1 19971 : [TsalUj 1198^ : 
Nakamura,.2000, : .Hiorth & Williams^2010,) , in that the coUisionless 



Boltzmann equation (CBE) is included as an additional constraint 
on the entropy functional. 

In this work, the aim is to study the lowest-order solutions of 
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Figure 1. The framework of the equihbrium statistical mechanics for self-gravitating systems. 



the statistical-mechanical theory we developed in iHel (l2012h . The 
paper is organized as follows. In Section |2] the framework of the 
statistical mechanical is briefly reviewed. Section |3] presents the 
solutions of the second-order equations for the equilibrium states. 
A summary and the conclusions of our study are given in Section|4] 



2 FRAMEWORK OF THE STATISTICAL MECHANICS 
AND ITS SECOND-ORDER EQUATIONS OF THE 
STATISTICAL EQUILIBRIUM STATE 

2.1 Framework of the statistical mechanics 

We outline the framework of the statistical mechanics of self- 
gr avitating s ystems in Fig.[T] The details of the theory can be found 



2.1.1 Local ergodicity and the Boltzinann-Gibbs entropy 

Ergodicity, or the equiprobability principle, is fundamental for 
equilibrium statistical mechanics. However, we demonstrated 
that the ergodic hypothesis is not valid for the whole self- 
gravitating system. It seems that ergodicity breaking may 
be a common feature for long-ran ge interaction systems 
(e.g. l Mukamel, Ruffo & Schrei ber 2005 ; Campa, Dauxois & Ruffd 
|2009|) . Fortunately we can re-establish the local ergodicity in any 
local finitely small volume element of the system by treating grav- 
itating particles as indistinguishable. With the validity of local er- 
godicity, we can still define entropy with the Boltzmann relation 
S = In Vl^, with W being the Fermi-Dirac-like microstate number 
as 



uj{x, v)\ 



n{x, v)\{uj{x, v) — n{x, v))! 



(1) 



where n{x,v) is the velocity distribution in the volume element, 
and Lu{x,v) — A^'xA'^v/h'^ denotes the degeneracy at the velocity 
level V for the macrocell, with hg being the volume of the minimum 



phase-space element (i.e. microcell). As a result, the entropy-based 
equilibrium statistical mechanics for self-gravitating systems can 
be established. 

The kinetic-theoretical equation for self-gravitating systems is 
the Boltzmann equation: 
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(2) 



Here, f{t,x,v) is Ihe fine-grained distribution function (DF), and 
r[/] is the encounter operator, given by 



r[f(wi,t)] = N / dVa 



9xi 



{wi,W2,t) 



(3) 



where w — (x, v) denotes the phase-space coordinates, g is the 
two-body correlation function, and $12 is the gravitational po- 
tential between the two bodies 1 and 2 (see equation 7.47 of 
iBinnev & Tremainell2008l) . 

If either 8^12/8x1 or g is vanishing, so that r[/] = 0, then 
the non-degenerate condition n{x,v) /uj{x,v) <^ 1 holds. With 
this non-degenerate condition, from equation ([T), we can derive the 
Boltzmann-Gibbs entropy: 



S[F] 



F(x,v)lnF{x,v)d^xd\. 



(4) 



where F{x,v) = n{x,v) /Ltj{x,v) is the coarse-grained DF. How- 
ever, if r[/] 7^ 0, in this case the non-degenerate condition is in- 
valid, and the coiTect entropy should be the Fermi-Dirac form as 



S[F] 



d^xd^v[F{x,v)lnF{x,v) 



+ (1 - F{x,v))\n{l - F{x,v))]. 



(5) 

Because the CBE (i.e. equation |2] with r[/] — 0, also 
the Vlasov equation) is a good approximation for self-gravitating 
systems if the particle number A'' of a system is very large, 
the Boltzmann-Gibbs entropy s hould be the proper entropy form 
for se lf-gravitating systems (cf. iTremaine. Henon & Lvnden-Belll 
[l98^ . 
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The CBE might not be a good approximation in very dense 
regions, such as the center of a simulated A'^-body self-gravitating 
dark matter halo, where two-body relaxation might play an impor- 
tant role. In this case, the CB E would be better replaced by, say , 
the Balescu-Lenard equation ( iHevvaertsI I20T0I ; |ch avani sl l201(]|) . 
and the Fermi-Dirac-type entropy of equation (|5), rather than the 
Boltzmann-Gibbs entropy, may be the correct entropy form. 



2.1.2 Additional macroscopic constraints 

The statistical mechanics of self-gravitating systems also requires 
that the CBE should be considered in some way as the addi- 
tiona l macroscopic constraint besides mass and energy conserva- 
tion jMadsen|[l987h , and'He (2012) has developed a way to prop- 
erly manipulate this additional constraint for spherically symmetric 
self-gravitating systems. First, we integrate the steady-state CBE 
(i.e. df/dt — in equation |2] with r[/] = 0) to obtain its v ari- 
ous o rders of velocity-moment equations (see equation All of lHel 
l20Tl) . 

dr r r (m + 1) 



,d<l> 



dr 



(6) 



with the even integers m,k ^ 0. These moment equations are then 
tr anslated in to their equivalent virialization forms (see equation 9 
of lH3l2012h as 



-2{m^ - l)Tk+2.m + 2(k + l)(m + 2)Tk,m+2 
+ (fc + l)(m + l)Wk,m = Mm + l)r^Pfc+2,. 



(7) 



These generalized virial equations are just the proper forms for 
the additional macroscopic constraints besides mass and energy 
conservation, which is an extension to the original prescription by 
[white & Narayan ( 1987). 

Finally, we emphasis that the coarse-grained DF F(x, v) and 
the fine-grained DF f{x, v) are different, in that the entropy func- 
tional, equation © or (|5), is constructed with F, but the DF in the 
CBE is /. There is no simple relationship between F and /, and 
hence we should not directly include the steady-state CBE as the 
constraint in place of the series of generalized virial equations. 



2.1.3 Out-of-equilibrium stationary states and the statistical 
equilibrium state 

As outlined above, the equilibrium states of self-gravitating sys- 
tems consist of both mechanical and statistical equilibria, the for- 
mer being characterized by a series of velocity-moment equa- 
tions {6]l, and the latter by the statistical equilibrium equations, 
which should be derived from the entropy principle, &S — 0, with 
all the relevant macroscopic constraints being considered. 

For a realistic gravitating system, however, such equilibrium 
states may not always be attained. Usually, if the collisionless re- 
laxations, for example violent relaxation, do not last long enough, 
then the s ystem will be trapped in so-called quasi-stationary states 
(see ICamp a et al. 2009; Chavanis .2006), whose lifetime depends 
on the total particle number A'^, proportional to N/ In TV. In the 
Vlasov limit, namely N ^ 00, this not-fully-relaxed system will 
finally be frozen into the out-of-equilibrium stable stationary state 
(i.e. SS 7^ 0) and the statistical equilibrium state is never reached. 

The statistical equilibrium state, characterized by the entropy 



SS ^ 



5S > 



Figure 2. Illustration of existence of both the stationary states and the sta- 
tistical equilibrium state. The thin lines indicate the stationary states, while 
the thick line indicates the statistical equilibrium state. The statistical equi- 
librium state is picked out by the entro py principle {5S = 0), which is 
exactly the saddle-point entropy state jHe & Kan J201(ll201 ih . whereas all 
the other stationary states do not satisfy the entropy principle {SS 7^ 0). 



principle SS = 0, is just a limit state, into which all the station- 
ary states will converge, assuming that the collisionless relaxation 
lasts long enough. We depict the relationship between the station- 
ary states and the statistical equilibrium state in Fig. |2] 

In reality, when a gravitating system, say a dark matter halo, 
experiences a monolit hic collapse from the cold initial conditions 
/Wana & White" 2009) or repeated mergers of progenitor haloes 
(Syer & White 1998), it will, according to our viewpoint, attain the 
statistical equilibrium state (or at least reside in a state very close 
to this equilibrium state), whose mass-den sity profiles are described 
by well-known empirical universal law s dNavarro. Frenk & White! 
'l997l: lMoore et al.ll999l ; lEinastjl96^ . If the systems are in quasi- 
stationary or stationary states, they may exhibit any p ossible den- 

Albadal 



sity shapes, which differ from the universal laws (e.g. I 
1985; 'Jing & Sutdl2OO0l ; Izhanell2004l ; lOin. Shan & Tilquiij|2008l; 

Lemze et al.l l20Ilh . 

Incidentally, the lLvnden-Belll jl967h approach turns out to be 
a good way to explain qualitati vely and quantitative ly some features 
of quasi-stationary states (see ICampa et al.ll2009l and references 
therein). Statistical theories that quantitatively describe these quasi- 
station a ry states have been const ructed (see Levin, Pakter & Teles 
'2008a|; iLevin. Pakter & Rizzatol boOSbl ; iRizzato, Pakter & Levii] 
i2009i) . 



2.2 Second-order equations of the equilibrium state 

We can express the coarse-grained DF F formally as the Taylor 
expansion 



F{X,V) = exp - ^ \k,m,n{x)viV^v'^ , 



with 



1 d 



fe-|-m + n 



lnF{x,v) 



k\mln\ dv'^dv^dv" 



(8) 



(9) 
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where \k.m,nix) are functions of the spatial coordinates x, and 
these undetermined coefficients should be determined by the 
joint mechanical and statistical equilibrium equations, and fc, m, n 
should be even integers, as there is no global momentum and angu- 
lar momentum for the system. 

We may expect that if all the expansion coefficients of F{x, v) 
in equation ([8} are exactly detemined, the predictions will give 
accurate agreement with the simulation results. For realistic opera- 
tions, however, we have to make some truncation to the DF. 

In this work, we truncate F{x, v) of the spherically symmetric 
system to the second order (i.e. the lowest order) and we directly 
replace the expansion coefficients by the corresponding moments 
as 



F{r,Vr,ve,v^) = 



(27r) 2 cr,.<T; 



Here, the density p and the velocity dispersions cr^ and are all 
functions of r. Because of the spherical symmetry, the dispersions 
of vg and are the same, so we uniformly use a'^ to represent 
oru2. 

With this truncated DF, we use our theory of the statistical me- 
chanics to derive the lowest-order (i.e. the second-order) equations 
as 



dpr- , 2 
dr r 
dlnp dlnpt 



-pt) 



dr dr 
3 d In p 1 dp,. 
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— 47rr 



2 _ ^Gm 



V? denotes both pg and 



(11) 



(12) 



(13) 



(14) 



with 



where pr = pv^, and pt = pv^ 
A in equation l ll3b being t he Lagrangian multiplier. Equation dllt 
is just the Jeans equation teinnev & Tremaine 2008), which is the 
lowest-order moment equation of equation (|6j (i.e. m — k = 0). 
Equations il2\ and ( 113b are two e ntropy stat ionary equations, de- 
rived from the entropy principle in lHel(l2012h . Equation (1141 is the 
differential form of the mass function, m(r) = Anr^pdr. See 
section 4.3 of lHel j2012l) for details about these equations. 



3 RESULTS 

The purpose of developing the statistical mechanics of self- 
gravitating systems is to provide a physical account of all the em- 
pirical relationships found in cosmological simulations concern- 
ing dark matter haloes. The above equations, however, are just 
the lowest-order approximations, derived from the truncated DF 
of equation ([8}. It should not be expected that these lowest-order 
equations will provide a good match to the simulation data. Never- 
theless, these lowest approximations may still reveal some features 
of this statistical mechanics, and also be valuable in contributing to 
an understanding of the structure of self-gravitating dark haloes. 

With a dimensional analysis, the Lagrangian multiplier A in 
equation il3i can be written as 



A = 



(15) 



inGpsrl ' 

where G is the Newtonian constant and ^ is a dimensionless num- 
ber. In all reasonable solutions, A should always be negative, and 
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Figure 3. Density profiles of the isothermal solutions of equations i ll It - 
(14). The behavior of the solutions depends on the initial values of pi in 
the inner regions, but all the solutions converge to the singular isothermal 
sphere p oc r^^ in the outskirts of the system. 



hence p is positive. We can always set = 1 by choosing appropri- 
ate values for the characteristic density ps and scale r^. The initial 
radius is chosen as rt = 10~'^rs. Given the initial values, we nu- 
merically solve these equations with the fourth-order Runge-Kutta 
method. 

According to the initial conditions, we can classify the solu- 
tions into the following three types. 



3.1 Isothermal solutions 



We define an effective velocity dispersion ctq as 



(Jq = 2nGpsr^- 



(16) 



So if = at — (Tq, the solutions of equations illi - il4t are 
isothermal solutions, but the solutions of density p(r) also depend 
on its initial values. For instance, if the initial value of pi at the 
radius r^ is set as /2-KGr1 , the equations will give the solution 
of the singular isothermal sphere p oc r~^, while if pi is set by the 
lEinasto ( 1965) profile 

ln(p(r)/p,) = --[(r/r.)"-l], (17) 

Q 

then the solution is flattened in the inner region and then converges 
to a singular isothermal sphere in the outskirts. We present the 

above two solutions in Fig. |3] 

Such a solution, similar to that of iLvnden-Belll ( Il967h . has 
infinite mass, energy and spatial extent, and hence is not realis- 
tic. Moreover, it is well known that such an isothermal solution is 
thermally unstable, which will lead to gra v othermal collapse (or 
gravothermal catastrophe) (lAntonov|[T962t iLvnden-Bell & Wood 
119681) . 



3.2 Divergent solutions 

We define the ratios of the radial and tangential velocity dispersion 
over the effective velocity dispersion (Tq as 



Ct 



2 



(18) 
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Figure 4. The divergent solutions of equations dill - )14t . The lEinastel j 196^ density profile is used to set the initial value of pi. The left two panels are for 
the case of Cr = 1-2 and Ct = 1, and the right panels are for Cr = 1 and Ct = 1.2. The dashed lines in panels (a) and (c) indicate the singular isothermal 
sphere, p oc r~^. 




Figure 5. The convergent solutions of equations 11 U - (14). The lEinastol i 19651) density profile is used to set the initial value of pi . The left two panels are for 
the case of Cr < 1, and the right panels are for Ct < 1. The dotted lines at the ends of the density functions indicate the truncation caused by the vanishing 
radial velocity dispersion . Different initial values of Cr and Ct are indicated in the figure. The dashed lines in panels (a) and (c) indicate the singular 
isothermal sphere, p oc r~^. 
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So for the isothermal solution, Cr = Ct = 1. Now we treat the 
initial Cr and Ct at as two free parameters, and see how the dif- 
ferent initial values of Cr and Ct will affect the behavior of the 
solutions. In Fig. [4] we present two solutions with either Cr or Ct 
being larger than 1. We can see that the density functions of both 
the solutions go outwards with respect to the radius r at a shallower 
slope than that of the singular isothermal sphere, and that all the ve- 
locity dispersions increase with r. So these solutions are divergent 
and unconfined, and of course are even more thermally unstable 
than the isothermal solutions. 



3.3 Convergent solutions 

If either Cr or Ct is smaller than 1, the solutions will be convergent 
and finite. In Fig. [5] we show these solutions with different initial 
values of Cr and Ct. The density functions at the outer bound- 
ary of the system are automatically truncated as a result of the 
vanishing radial velocity disper sions. Although the results cannot 
match the simulations well (e.g. lNavarro et al.]|2010l : lLudlow et al.l 
these lowest-order solutions do possess the following prop- 
erties: (i) the density profile is finite, which does not lead to a state 
with infinite mass or energy; (ii) the velocity dispersions are vari- 
able functions of the radius r, and (iii) the velocity distributions 
are anisotropic in different directions. The last two points indicate 
that the statistical equilibrium of self-gravitating systems is by no 
means the thermodynamic equilibrium. 

We use the Boltzmann-Gibbs entropy in our statistical- 
mechanical theory. It might be a misunderstanding that the 
Boltzmann-Gibbs statistics must lead to the Maxwell-Boltzmann 
velocity distribution. However, by including the additional con- 
straints of the generalized virial equations (|7), the Boltzmann- 
Gibbs statistics is endowed with new properties, in that the theory 
may give any possible velocity distributions. 

As addressed previously, these convergent solutions are just 
the lowest-order approximations, derived from the truncations of 
the DF, equation ([8]l, and its moment equations As a result, 
these lowest-order equations cannot provide a good match to the 
simulation data. We are optimistic in expecting that higher-order 
solutions of the statistical-mechanical theory will give much better 
agreement with the realistic data. 



tionary states will converge, if that the collisionless relaxation lasts 
long enough. 

The coarse-grained DF is formally expressed as a Taylor se- 
ries, and if all the expansion coefficients are determined, the the- 
ory will give exact predictions. For practical purposes, however, we 
have to truncate this Taylor series. In this work, we truncate the DF 
for the spherically symmetric system to the second order (i.e. the 
lowest order). With this truncated DF, we use our theory to derive 
a set of second-order equations, and it was our aim to solve these 
equations in this paper. 

According to the initial conditions, we can classify the solu- 
tions into isothermal, divergent and converge nt solutions. We fin d 
that the isothermal solutions, similar to that o f lLvnden-Belll ( ll967h . 
have infinite mass, energy and spatial extent, and hence are not ac- 
ceptable. Moreover, it is well known that such isothermal solutions 
are thermally unstable, which will lead to gravothermal collapse for 
the system. The divergent solutions should be even more thermally 
unstable than the isothermal solution. 

The convergent solutions seem, however, to be reasonable. Al- 
though the results cannot match the simulations well, these lowest- 
order solutions show that the density profiles are finite and con- 
fined, the velocity dispersions are variable functions of r, and the 
velocity distributions are also anisotropic. In addition, the solutions 
indicate that the statistical equilibrium of self-gravitating systems 
is by no means the themodynamic equilibrium. 

Our statistical-mechanical theory is based on the Boltzmann- 
Gibbs entropy. By including the generalized virial equations as ad- 
ditional macroscopic constraints, this new Boltzmann-Gibbs statis- 
tics may give any possible velocity distributions. These convergent 
solutions are just the lowest-order approximations, derived from 
the truncations of the DF and its moment equations, but they have 
already manifested the qualitative success of our theory. We expect 
that higher-order solutions of the statistical-mechanical theory will 
give much better agreement with the simulation results for dark 
matter haloes. 
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4 SUMMARY AND CONCLUSIONS 

InlH3 j2012h . a systematic theoretical framework was formulated 
for the statistical mechanics of spherically symmetric collision- 
less self-gravitating systems based on the Boltzmann-Gibbs en- 
tropy. According to the statistical-mechanical theory, the equilib- 
rium states of self-gravitating systems consist of both mechanical 
and statistical equilibria. The mechanical equilibria are character- 
ized by the steady-state Vlasov equation, or its equivalent velocity- 
moment equations, and the statistical equilibrium equations should 
be derived from the entropy principle. 

For a realistic gravitating system, however, such equilibrium 
states may not always be attained. Usually, if the collisionless re- 
laxations do not last long enough, the system will be trapped in 
quasi-stationary states, whose lifetimes depend on the total parti- 
cle number. In the Vlasov limit, this not-fully-relaxed system will 
finally be frozen into the out-of-equilibrium stable stationary state, 
and the statistical equilibrium state is never reached. Hence, the 
statistical equilibrium state is just a limit state, to which all the sta- 
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